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Abstract. The transport properties of the planar rotator model on a square lattice 
are analyzed by means of microcanonical and non-equilibrium simulations. Well 
below the Kosterlitz-Thouless-Berezinskii transition temperature, both approaches 
consistently indicate that the energy current autocorrelation displays a long-time 
tail decaying as t . This yields a thermal conductivity coefficient which diverges 
logarithmically with the lattice size. Conversely, conductivity is found to be finite in the 
high-temperature disordered phase. Simulations close to the transition temperature 
are insted limited by slow convergence that is presumably due to the slow kinetics of 
vortex pairs. 
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1. Introduction 

Physical phenomena in reduced spatial dimension (d = 1,2) are often qualitatively 
different from their three-dimensional counterparts. The overwhelming role of statistical 
fluctuations and the presence of constraints in the motion of excitations can lead to 
peculiar effects like the impossibility of long-range order. In the context of non- 
equilibrium statistical mechanics, the existence of long-time tails in fluids [I] leads to 
ill-defined transport coefficients, thus implying a breakdown of the phenomenological 
constitutive laws of hydrodynamics [2] . 

A remarkable example is the anomalous behavior of heat conductivity for d < 2. 
This issue attracted a renovated interest within the statistical-mechanics community 
after the discovery that the thermal conductivity of anharmonic chains diverge in the 
thermodynamic limit Since then, those anomalies have been detected in a series of 
different models. An exhaustive account is given in Ref. [3], where the effects of lattice 
dimensionality, disorder and external fields for the validity of Fourier's law are discussed 
in detail. The signature of anomalous behavior is a non-integrable algebraic decay of 
the correlator of the heat current J (the Green-Kubo integrand) at large times 

(J(t)-J(O)) oc r (M , t^+oo (1) 

where < a < 1 and ( ) is the equilibrium average. For a finite system of linear size 
L this implies that the finite-size conductivity k(L) diverges in the L — > oo limit. In 
fact, in the framework of linear- response theory, k can be estimated by cutting-off the 
integral in the Green-Kubo formula at the transit time L/v (v being some propagation 
velocity of energy carriers). Taking into account Eq. (0) one straightforwardly obtains 
k oc L a . 

Simulation studies of specific models |lj as well as analytic arguments lead to the 
surmise that anomalous conductivity should occur generically whenever momentum is 
conserved. Moreover, the exponent a should be largely independent on the microscopic 
details as suggested by a the renormalization-group calculation of Ref. [Sj that predicts 
a = (2 — d)/ (2 + d) . For d = 1 the resulting value a = 1/3 is roughly close to the 
numerical estimates although some substantial deviations have been observed in specific 
cases [3 El IE]- The situation is even more controversial d = 2 where the predicted 
decay yields a logarithmic singularity which is consistent with simulation data jH] , while 
other works report significative deviations and dimensional crossovers jl()j . 

When extending the analysis to the 2d case, one naturally wonders how the 
possibility of observing critical phases may affect the anomalous energy conduction. 
The first example one may think of is of course the Ising model. It has however been 
shown that, at least for a specific choice of the spin dynamics, the latter displays a 
normal conductivity at all temperatures [TT]. In a more general perspective, this is 
consistent with the idea the breakdown of momentum conservation removes transport 
anomalies. Indeed, the field-theoretic counterpart of the Ising model, the so called 4 - 
theory, acquires an on-site non-linear interaction that breaks translational invar iance. 
Although we are not aware of any study of this model in the ordered phase, it is known 
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that the lattice 4 -model in the high-temperature phase displays a finite conductivity 
in the thermodynamic limit [T2"] . 

In this present paper, we present a simulation study of the transport properties of 
a model of rotators coupled on a square lattice, akin to the celebrated XY-model (see 
e.g. PS] and references therein for a comprehensive review). As it is well known, the 
latter is characterized by the presence of the so called Kosterlitz-Thouless-Berezinskii 
(KTB) phase transition at finite temperature between a disordered high-temperature 
phase and a low-temperature one, where vortexes condensate. 

As recalled above, the fact the momentum (actually the angular momentum) is a 
constant of motion, makes this model a candidate for observing anomalous behavior. 
On the other hand, its Id version is the only known exception to this requirement and 
displays normal trasport, due to the presence of "dynamical defects" in the form of 
localized rotations that act as scattering centers for the heat carriers |14j . On the basis 
of this observation it is extremely interesting to investigate if and how the vortices play 
a similar role in the 2d case. 

Before entering the details of the present work, it is important to mention that 
some evidence of the role of the vortex unbinding on the transport properties of the 
(modified) XY-model have been reported in Ref. j!5j . Nonetheless, those results are 
mainly qualitative and we are thus motivated to undertake a more detailed analysis. 

The paper is organized as follows. In section 2 we introduce the model and its 
microcanonical simulation. In Section 3 we recall the technique we used to investigate 
the non-equilibrium stationary state. The outcomes of numerical simulations for the 
disordered and critical phases are reported in Sections 4 and 5, respectively. Finally, we 
summarize and discuss our results in Section 6. 

2. Hamiltonian dynamics of the XY model 

The XY or planar rotator model consists of a set of classical "spins" S r of unit length 
confined in a plane, whose orientation is specified by the angle 9 r , with r = (i, j) being 
an integer vector labelling the sites of a square lattice of size N = N x x N y . It is 
known that this model does not admit equations of motion and therefore its canonical 
dynamics is usually simulated either by Monte-Carlo methods [To! ITTj or by Langevin 
type equations [THIEI. Microcanonical approaches consist instead, either in considering 
a three component spin model |2P. or into adding a kinetic energy term [21] . The latter 
method, which we follow in the present work, can be also generalized to other physical 
systems (see e.g. the application to lattice gauge theories [22]). All these different 
dynamics should display the same static properties, as it has been verified up to some 
extent jSUESj- We thus consider the Hamiltonian 

^ = E#+E[i- «*(^ " Or)) , (2) 

r Z (r,r'> 
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where p r = 8 r is the angular momentum of the rotator. The sum ranges over the four 
nearest neighbors of site r, namely r' = r ± x and r' = r ± y where x and y are the unit 
vectors parallel to the lattice axis. It could be shown that (j2J) is obtained as the classical 
limit of a quantum Heisenberg Hamiltonian with an anisotropy term ^ r (S r 2 ) 2 , using the 
representation introduced in Ref. [23]. We have set both the inertia of the rotators and 
the ferromagnetic coupling constant to unity so that the only physical control parameter 
is the energy per spin e = TC/N. Actually, there is a second constant of the motion, the 
total angular momentum P = J2rPr, whose choice affects the results in a trivial way. In 
the numerical simulations we set P = to avoid global ballistic rotation. 

A previous study of the static properties of (J2J) j2H] showed that the system 
undergoes a KTB transition [T3] at e = &ktb ~ 1-0, corresponding to a kinetic 
temperature Tktb ~ 0.89, which is in agreement with the value 0.894(5) obtained 
in the canonical ensemble flj. One of the most striking features of the XY model is the 
presence of strong finite-size effects j2H], e.g. the existence of a sizable magnetization 
for large samples, despite the fact that long-range order cannot occur for the infinite 
system. This has also some consequence on the dynamical correlation of the finite-size 
magnetization j2E]- 

In the framework of linear-response theory, heat transport properties can be 
analyzed by computing the autocorrelation function (or, equivalently, the power 
spectrum) of the total heat current vector J at equilibrium. We thus need a microscopic 
expression that can be worked out by the procedure followed for other similar models 
]3]. In brief, it amounts to writing down a discretized continuity equation and, by means 
of the equation of motion, identify the proper expression of the local flux in terms of 
the canonical variables (9 r ,p r ). For model (J2J), J = (J x ,J y ) can be written as a sum 
over all lattice sites 

JX = \ E Sin(0 r+ x " 0r) pr+x + Or] (3) 
^ r 

p = ~£sin(0 r+ *-0 r ) \e r+ y + e r ] (4) 



r 



This latter expression is the correct one in the microcanonical ensemble with P = 0. 
Incidentally, notice that a suitable counterterm should be subtracted out if one wishes 
to work in a different statistical ensemble [27] . 

The numerical integration of the equations of motion (with periodic boundary 
conditions) is performed using the fourth-order McLahlan-Atela algorithm j2B] , which is 
an explicit scheme constructed from a suitable truncation of the evolution operator that 
preserves the Hamiltonian structure. One of the major merits of symplectic algorithms 
is that the error on the energy does not increase with the length of the run. The chosen 
time step (0.01-0.05 in our units) ensures that in every simulation energy fluctuates 
around the prescribed value with a relative accuracy below 10~ 5 . 

As mentioned above, the main quantity of interest is the flux autocorrelation 
function (J(t) ■ J(0)). For numerical purposes, we find more convenient to evaluate the 
power spectra S(f), i.e. the squared modulus of the Fourier transform of each component 
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of J, averaged over a set of different initial conditions. These initial conditions were 
chosen by letting 6 r = and drawing the 6 r = at random from a Gaussian distribution 
with zero average and unit variance. The momenta are then all rescaled by a suitable 
factor to yield the desired value of the total energy. A transient is elapsed in order to 
start the averaging from a more generic phase-space point. 



3. Non— equilibrium simulations 

The non-equilibrium simulations have been performed by coupling all the rotators on 
the left and right edges of the lattice with two thermal baths operating at different 
temperatures T + and T_. Periodic and fixed boundary conditions have been adopted 
in the direction perpendicular (y) and parallel (x) to the thermal gradient, respectively. 
Thermal baths have been simulated by applying the Nose-Hoover method: 

dV 

Or = - W -0r[C J + S l , 1 + C J S hN } (5) 







0^ 1 


\k B T + 






eT 1 


\k B T- 



- 1 



Here, V is the potential associated with (J2J), 0± are the thermostats' response times, 
and 5 is the usual Kronecker symbol. Notice that each rotator is thermostatted 
independently and, accordingly, the Nose-Hoover variables are vectors of length 
N y . For computational purposes, simulations have been performed by fixing the aspect 
ratio R = N y /N x < 1. The choice of R results from a trade-off between minimizing the 
number of rotators (R small) and dealing with a genuinely 2d lattice (R ~ 1). Indeed, 
too small values of R would require considering larger system sizes to clearly observe 2c? 
features. For small lattices, we checked that the results are almost independent of the 
ratios employed hereby. 

Since we are interested in the average values of the flux we have to check that the 
non-equilibrium stationary state is indeed attained. To this aim, we monitored that the 
average fluxes towards the boundaries 

J + =-Y,0h (6) 



j 

52 



j 

were equal to the flux in the bulk, namely J ± = J x (the overline denotes a time average 
henceforth). 

As observed before jlj, the choice of the thermostat response times is crucial to 
the time needed to reach the stationary state, and to control the values of thermal 
resistance at the boundaries. In order to fasten the convergence, the initial conditions 
have been generated by thermostatting each particle to yield a linear temperature profile 
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along the x direction. This method is very efficient, especially for large lattices, when 
thermalization within the bulk may be significantly slow. 

Once the steady non-equilibrium state is attained, the relevant observables are 
accumulated and averaged in time. In particular we evaluated the local kinetic 



formulae (j3J), as well as the average fluxes towards the reservoirs J ± . For obvious 
simmetry reasons, we expect (and indeed found) the J y vanishes up to the statistical 
accuracy. Moreover, T r depends only on the location along the x axis, and a further 
average is performed along y. Once the the average flux is computed, we evaluated the 
thermal conductivity coefficient from Fourier's law as 



The last equality is only approximate, since the actual thermal gradient within the 
lattice is usually smaller than (T+ — T-)/N x , due to boundary resistance effects jlj. In 
other words, the coefficient evaluated in this way should be regarded as an effective 
conductivity including both boundary and bulk scattering. As we are going to show 
in the next section the rescaled stationary temperature profiles, obtained for different 
values of N x are such that the temperature gradient actually scales like iV" 1 . Therefore, 
definition (JZj) yields the same scaling behavior of the bulk conductivity. 

4. The disordered phase 

Let us start discussing the case of high energies or temperatures where the system 
is away from criticality. In the non-equilibrium simulations we fixed T + = 1.5 and 
T_ = 1.4 which are both well above the KTB transtion temperature. Both response 
times of the Nose-Hoover thermostats have been set to the same value, ©± = 1.0. 
In fact, we have found empirically that such a choice minimizes boundary impedance 
effects, thus allowing for larger values of the heat flux. In order to improve the statistics, 
the temperature profiles and the measures of the heat flux have been averaged over 16 
independent initial conditions. Numerical simulations have been performed by fixing 
the value of the apsect-ratio to R = 1/2 and by increasing N x up to 140. Some of 
the temperature profiles are reported in Fig. They all exhibit a linear shape, which 
testifies to the expected temperature profile when Fourier's law holds. In Fig.|2]we show 
that the thermal conductivity, as defined by (J7J), is independent of N x , with fluctuations 
around the average value, extending up to some 10 %. 

These results have been compared with those obtained from equilibrium 
simulations. The value of the energy density has been chosen e = 2. Actually, this 
value corresponds approximately to the average temperature, (T + +T_)/2, of the above 
mentioned non-equilibrium simulations. We want to point out that, for very large values 
of e, the kinetic energy dominates over the potential one and the system approaches 
the integrable limit of independent free rotators. Accordingly, in this limit the lattice 
is expected to behave as a perfect insulator, since the time scale for transmitting any 



temperature hsT? 



6% and the average energy fluxes J x and J y , as defined by 




(7) 
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Figure 1. The temperature profile in the high-temperature phase for different values 
of N x , which have been rescaled to the unit length. The response times of the Nose- 
Hoover thermostats are 0± = 1.0, which guarantee negligible boundary impedance 
effects. 
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Figure 2. Finite-size conductivity in the disordered phase. Nose-Hoover thermostat 
with response times fixed to 1.0; Each point results from an average of 16 independent 
runs of about 10 6 time units each. The error bars are the error on the mean and the 
horizontal line is the average of all the measured values, k = 0.24. 



energy fluctuation diverges. In this respect, the choice e = 2 is appropriate, also because 
one can observe convergence of the quantities of interest over reasonable simulation 
times (typically, 10 6 time units). In Fig. 01 we show the heat flux spectra: they are 
independent of the lattice sizes and tend to a constant for small frequencies. This in 
a clear confirmation that no long-time tail is detectable and the thermal conductivity 
is a well-defined quantit in the thermodynamic limit. It should be remarked that the 
lineshape of these spectra cannot be fitted by a simple Lorenzian. This implies that in 
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the high-temperture phase the decay of time-correlations cannot be reduced to a simple 
exponential law. 




f 



Figure 3. Power spectra of the heat current in the disordered phase for three different 
lattice sizes, N = 8 2 , 16 2 ,32 2 . Data are averaged over 200 random initial conditions. 
In order to minimize statistical fluctuations, a further averaging of the data over 
contiguous frequency intervals has been performed. The three spectra actually almost 
overlap: this is why we have presented them after a vertical arbitary shift in order to 
better distinguish one from each other. 



5. The critical phase 

A more interesting situation appears in the low-temperature phase. For what concerns 
the non-equilibrium calculations, we have fixed the thermostats' temperatures T± to 
be well below the KTB transition value. At variance with the high-temperature 
phase, here the values of the response time of the Nose-Hoover thermostats have to 
be properly tuned in order to minimize boundary impedence effects. In particular, we 
have determined empirically the values + = 2.0 and G_ = 6.0. Such different values, 
between themselves and also with respect to the high-temperature phase, indicate that 
fluctuations have to be slowed-down significantly in order to cope with the typical time 
scales of the dynamics. Moreover, we have performed the same statistical averaging as 
in the high-temperature case. 

In Fig. H we show the temperature profiles obtained for T + = 0.5, T_ = 0.4 and 
different values of N x . They exhibit a good data collapse when N x is rescaled to the 
unit length. This confirms that also in the low-temperature phase the thermal gradient 
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scales like the inverse of the system size, VT ~ (T + — TJ)/N X . On the other hand, 
the temperature profile has assumed the typical non-linear shape, which testifies to 
anomalous thermal conductivity. Moreover, this shape is similar to the temperature 
profiles of the one and two-dimensional Fermi-Past a-Ulam model jl]. 

0.5 



0.48 
0.46 
0.44 
0.42 

°- 4 0.2 0.4 0.6 0.8 1 

i/N 

x 

Figure 4. The temperature profile in the low-temperature phase for different values 
of N x , which have been rescaled to the unit length. The response times of the Nose- 
Hoover thermostats are G+ = 2.0 and 0_ = 6.0, which minimize boundary impedence 
effects. 

The finite-size thermal conductivity as a function of the longitudinal size N x 
is reported in Fig. |3] for three different values of the boundary temperatures (with 
T + — T_ kept fixed to 0.1). As expected, increasing T± the conductivity decrease. More 
importantly, for fixed temperatures, the data all exhibit a systematic increase with N x . 
In analogy with what found in the 2d Fermi-Pasta-Ulam model the data for lower 
temperatures (curves (a) and (b)) can be well fitted by a logarithmic law 

k(N x ) = d + C 2 \nN x . (8) 

Actually, a closer inspection of data set (c) reveals that the logarithmic fit is rather poor. 
In view also of the limited size range we have been able to explore, a convincing estimate 
of the growth law is unfeasible with the data at hand. This is presumably due to the fact 
that approaching Tktb requires longer times and sizes. In fact, we have observed that 
convergence of the averages considerably slows down in this region (about a factor 5 in 
passing from simulation (b) to (c)). The effect may be caused by slow thermalization of 
vortex pairs. Indeed, in this temperature region the vorticity starts to become sizeable 
[23] and a slower kinetics as well as relevant correction to scaling may thus be expected. 

Following the analysis performed for the high-temperature phase, we have 
performed also microcanonical simulations, with periodic boundary conditions imposed 
in both lattice directions. We only considered the energy density e = 0.5 which roughly 
corresponds to T = 0.45. 
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Figure 5. The thermal conductivity K versus N x in the low-temperature phase for 
T_ = 0.3, T+ = 0.4 (a), T_ = 0.4, T + = 0.5 (b) and T_ = 0.6, T+ = 0.7 (c). Each point 
is the result of an average over 16 independent initial conditions, each one lasting for 
10 6 time units. The error bars are of the order of the symbol's size. The solid lines 
are a best-fit with a logarithmic law, Eq. (JSJ. 

In Fig. |Hl we report the power spectra of the heat flux for different lattice sizes. 
For large frequencies (/ > f c ~ 10~ 3 ) we have a f~ 2 behaviour which suggests a fast 
(exponential) decay of the correlation at short times. The crossover frequency f c should 
be related to some typical time-scale for the hydrodynamic effects to set in. In the 
low-frequency limit, / < f c , the data are consistent with a logarithmic singularity of 
the type (see the inset of Fig. EJ) 

S(f) = A + B\nf (9) 

which, in turn, corresponds a t~ l tail of the autocorrelation function. According to 
the argument exposed below Eq. ([TJ, this would yield the logarithmic divergence (see 
Eq.fjHJ)). Altogether, we conclude that the equilibrium and nonequilibrium aproaches 
yield the same divergent behaviour. 

It must be admitted that the fitting of the low-frequency part with formula (jOJ) 
is convincing only for the larger sizes. Actually, a power-law fit S(f) oc f~ 0A is 
also compatible with the data obtained for smaller iV values. On the other hand, 
this estimate cannot be taken seriously for a twofold reason. First of all, it may be 
easily attributed to finite-size effects. Moreover, it would imply a power-law divergent 
conductivity which, in turn, would be inconsistent with the non-equilibrium data (see 
again curve (b) in Fig. EJ). 

In view of the above fact, one may wonder why equilibrium simulations should 
be much more sensitive to finite-size effects than nonequilibrium ones. A reasonable 
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Figure 6. Power spectra of the heat current in the low-temperature phase (e = 0.5). 
for four different lattice sizes N = 16 2 , 32 2 , 64 2 , 120 2 . Data are averaged over 400 
random initial conditions. In order to minimize statistical fluctuations, a further 
averaging of the data over contiguous frequency intervals has been performed. The 
spectra are presented after a vertical arbitrary shift to better distinguish one from 
each other. The inset is an enlargement of the low-frequency region in log-lin scale. 

qualitative explanation goes as follows. At low energies, where the isolated system 
is almost harmonic, it is very unlikely that a rotator turns from small oscillations to 
fast rotations. Moreover, for a small size system this process is even more unlikely 
as it demands a sufficiently large local energy fluctuation. Since this is the main 
scattering mechanism, one should wait for very long simulations before recovering the 
true asymptotic regime of energy transport. Conversely, in nonequilibrium simulations 
thermal baths act as external sources of fluctuations. These may favour the creation of 
the nonlinear excitations, thus shortening the time scale needed for the effects of the 
scattering process to be appreciated. 

6. Concluding remarks 

We have found numerical evidence that transport properties of the XY model on a finite 
lattice are drastically different in the high-temperature and in the low-temperature 
phases. In particular, thermal conductivity is finite in the former case, while in the 
latter it does not converge up to lattice sizes of order 10 4 . In the region where vorticity 
is negligible (T < 0.5) the available data suggest a logarithmic divergence with the 
system size analogous to the one observed for coupled oscillators [Oj. Close to Tkbt, 
where a sizeable density of bounded vortex pairs are thermally excited, our data still 
suggest a divergence, whose law we cannot reliaby estimate. 

We want to point out that these results have been obtained consistently for 
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equilibrium and nonequilibrium simulations. The equivalence between these approaches 
is not granted a priori. On the other hand, this property has been verified for many 
other similar models jlj, whose dynamics is of Hamiltonian type. In this respect, 
different choices of the dynamics (e.g. Monte-Carlo) may not necessarily lead to the 
same conclusion. 

This is a very important and interesting physical result: it indicates that in the 
low-temperature phase some materials or states of matter, may behave as anomalously 
efficient heat conductors. For instance, as a direct consequence of the studies performed 
in this paper on the 2d XY model, liquid Helium films should be included in this class of 
materials. An experimental test confirming the prediction of the logarithmic divergence 
of the heat conductivity with the system size would be highly appropriate and welcome. 

A complete hydrodynamic theory of the 2d XY model could certainly help in 
clarifying many of the aspects that our numerical approach cannot fully assess. For 
instance, an approach based on spin-waves proved out to be be effective in evaluating the 
dynamical correlations of the low-temperature phase [22] • On the other hand, estimates 
of the energy-current correlators may be technically more difficult. Indeed, the heat 
flux is a constant to leading order and, accordingly, one has to account for higher-order 
terms for the theory to make any sense. In this respect the calculations should be 
conceptually equivalent to estimating spin-waves lifetimes (HHl- It is not however clear 
how to include the effects of vortices and finite-size magnetization in this framework. 
The hydrodynamics in the high-temperature phase presumably should go through less 
technical troubles, although the construction of, say, a large-deviation functional is far 
from trivial. Anyway, an effort in this direction is in our future agenda. 
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